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Abstract 

We report on a study of the autocorrelation times of the local version of the Hybrid 
Monte Carlo (LHMC) algorithm for pure gauge SU(3). We compare LHMC to standard 
multi-hit Metropolis and to the global version of the same HMC. For every algorithm 
we measure the autocorrelation time (r) for a variety of observables and the string 
tension (cr) as a function of j3. The measurements performed on 8 and 16 4 lattices 
indicate that the autocorrelation time of LHMC is significantly shorter than for the 
other two algorithms. 
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Introduction 

In order to extract a meaningful number from a Monte Carlo (MC) simulation we have 
to associate a realistic error estimate to the measurement of the mean value. The ability of 
obtaining a good value for the error depends on our knowledge of the number of statistical 
independent samples that we have, which in turn implies a knowledge of the autocorrelation 
time r associated to the chosen algorithm. 

In lattice QCD, as we approach the continuum limit ( a — > 0, or (3 — > oo) we observe the 
well known phenomenon of critical slowing down. This means that, as a —>■ 0, r tends to 
diverge, thus requiring more and more iterations of a given algorithm to obtain an indepen- 
dent configuration. The approach to the continuum limit implies that the autocorrelation 
length of the system is diverging, or the lowest mass, expressed in lattice units is diverging. 
For a given lattice, a numerical simulation will be meaningful only if we limit ourselves to 
values of the coupling constant corresponding to correlation length significantly shorter than 
the lattice size, so for a given lattice the critical slowing down can be described by a critical 
exponent z parameterizing the dependence of r on the typical mass scale, r = amT z . 

In chapter 1 we give a brief overview of the HMC algorithm and the LHMC version, 
chapter two contains a description of the measurement of the autocorrelation time, chapter 
three describes in details the measurement performed, including a discussion on the choice 
of the optimal trajectory length. Results are presented and commented in chapter four. 

1 Hybrid Monte Carlo 

The HMC algorithm is one of a variety of Monte Carlo techniques to simulate statistical 
systems described by some local action 

We will concentrate on lattice gauge theories, therefore the action we are talking about 
is the usual Wilson's action, 

S p = ]T P[l - ^ReTr(U^n)U v (n + p)U^{n + + i/)E/_„(n + u))}. (1) 

X,fl,l> " 



The HMC prescription dictates to introduce fictitious conjugate momenta ir for every 
degree of freedom, and simulate the pseudo-Hamiltonian 



H(n,U) = Tr^(x) + S p (U). (2) 
The simulation begins by generating the momenta 7r's distributed according to the Gaussian: 

P(n) = exp('£T™l(x)), (3) 

then, starting from the randomly chosen momenta and the current C/'s a new set of 7r's 
and £7's are chosen as the end point of the trajectory generated by integrating Hamilton's 
equation of motion associated to the Hamiltonian (2). 

At the end of the trajectory, the small variation in energy introduced by the finite step 
size is corrected by a Metropolis-like step. 

All the details on this algorithm can be found in |1J B. The key point to observe is that, 
since all the degrees of freedom are evolved simultaneously, in order to preserve detailed 
balance, the Metropolis step must be applied to the whole lattice. This step consists in 
comparing the difference between the initial and the final value of the energy 

SH = H new — H i d (4) 

and the probability of accepting the new configuration is proportional to 

P « e~ SH . (5) 

8H is an extensive quantity and it will grow with some power of the volume. In turn, this 
implies that, in order to keep the acceptance constant, the step size will have to be reduced 
according to the same volume's law. 

1.1 Local Hybrid Monte Carlo 

Let us denote with Ue the gauge field variables sitting at even sites of the lattice and Uo those 
sitting at odd sites. The stochastic matrix describing the updating from the old variables 



(Ue, Uo) to the new variables {Ue, Uo) will be denoted by: 

W(U^,Th;U E ,Uo). (6) 

If we choose to update even sites first, then odd sites, this stochastic matrix factorizes in 

W(U E ', Uo; Ue, Uo)W(U^, Uo; U e , U ). (7) 

Since the Wilson action only involves nearest neighbor interactions, even sites are totally 
independent, so are odd ones and, as a result, each one of the factors of the above equation 
factorizes. For instance: 

W(Ue, Uo; U e , U q ) = \{ W{TT X , U Q ; U x , U ). (8) 

xeE 

Parallelism can be achieved easily by applying simultaneously all the elemental stochastic 
matrices to one checkerboard at a time. For each variable then, we apply the HMC scheme 
and the final Metropolis acceptance-rejection test only involves a local modification of the 
action therefore there is not any degrading of the acceptance due to the volume and we can 
use very large step sizes. 

The integration of Hamilton's equations of motion involves the computation of the force 
term which is made up by the staple of the link under investigation. This staple is constructed 
solely by links belonging to the opposite checkerboard, so the algorithm is equivalent to the 
integration of the equation of motion for a variable in an external field. This is, of course, an 
other clear advantage over the conventional HMC. In the global version of HMC, at every 
step the staple has to be computed anew, therefore the computational cost of a trajectory 
is a multiple of the cost of the first step. In the local version only the first step is expensive, 
successive steps are only a fraction of the cost that can be estimated around 10%. 

As we will see in the section devoted to the discussion of the results, the combination 
of these effects account for an inexpensive algorithm that moves effectively through phase 
space. 



2 Autocorrelation times 



For a set of n measures of an observable A we can define the autocorrelation function: 

n-\t\ 

C(t) =NJ2< {M- < A >){A i+t - <A>)> (9) 
i=i 

where < A > is the average over all the measurements and N is a constant fixed by the 
condition C(0) = 1. If we assume that, at large times t, this function decays exponentially 
C(t) ~ exp[—t/T exp \, then r exp provides us with a definition of the autocorrelation time. This 
way to compute r however can be affected by large errors, when r is small. 

This defines, of course, only the autocorrelation of the given observables. Different ones 
can in principle, and often do, have different autocorrelation times. To define the autocor- 
relation of a given algorithm, one should survey a complete set of observables and chose the 
highest r. This, being highly non practical, is approximated by a choice of selected operators 
and that will, operatively, define the r of the algorithm. 

Another definition of r is given by ||: 

1 oo 

Tint = 5 E C(t) (10) 
Z t=-oo 

However, this estimator of Ti nt present a whole new set of troubles because when the sum 
over t goes to \t\ » r the function C(t) contains noise much higher than the signal. Adding 
these terms we would over estimate the errors. The practical solution is a cut off in the sum 
in the previous equation, introducing a suitably constant M: 

i M 

Tint = -= E C(t) (11) 
Z t=-M 

In this way our calculation will be biased, since we are systematically neglecting terms in 
the sum (6). On the other hand we can systematically exploit that bias defining the sum 
extended over t's as long as C(t) > Our choice is q — 10. 

Whatever definition of r we choose, we noticed consistently that in order to have a 
reliable estimate, the length of our run must be approximately 1000 times the estimated 
autocorrelation time. 



Our r have been measured on a large set of Wilson loop, from 1 x 1 to 3 x 3. For these 
observables the behavior is not noticeably different, therefore the data shown in the following 
refer to the 3x3 Wilson loops. 

3 Tuning of the parameters 

There are two parameters characterizing HMC like algorithms, N s and St. It is well known 
that for the global HMC tuning of these parameters is crucial to improve on the autocor- 
relation of the algorithm. For the LHMC we have noticed that the autocorrelation time, 
of the monitored observables, depends very lightly on the parameters chosen. A value of 
T = N s 5t = 1 is nearly optimal in every case and all data shown have been obtained with T 
between 0.9 and 1.2, with N s varying from 2 to 4. These values correspond to an acceptance 
rate around 90%. 

When discussing autocorrelation times, these should be correctly given in CPU time and 
not sweeps but for the LHMC, as explained above, extra steps have only a marginal cost, so 
we give, for sake of simplicity, all the r's expressed in number of iterations. The same holds 
for multi-hit Metropolis, while in the case of global HMC, where different trajectories length 
have a very different cost, we limited ourselves to the fixed trajectory length of T = 0.3, with 
N s = 10 ||. All the simulat ions for multi-hit Metropolis were done with 8 hits per step. 

4 Results 

Autocorrelation times have been measured, besides for LHMC, for multi-hit Metropolis and 
the global HMC. In table 1 me show the autocorrelation times for the 3x3 Wilson loop 
on a 16 4 lattice, for all the three algorithms considered, while tables 2 and 3 concentrate on 
LHMC for different observables and /3's. 

From the data it is obvious that LHMC decorrelates much faster when measured in 
terms of sweeps, and even better if measured in terms of computational cost. This concept 
is emphasized in table 4 where we plot the time it takes to perform a number of sweeps 
corresponding to one r. For instance, at (3 = 6.0, if we assign, in arbitrary units, 1 to LHMC 



it takes 7.4 for multi-hit Metropolis and 400 for global HMC. 

Our goal was the measure of the dynamical critical exponent of the algorithm in the 
hope of confirming certain theoretical results obtained on the Gaussian model, claiming it 
to be one. This was not possible since we noticed an unexpected behavior of the r's. For 
every observables that we analyzed, we observed that, plotting r versus the autocorrelation 
length (or for that matter versus (3), we did not get a monotonically increasing function, as 
expected, but a peak at a finite value of the coupling. The position of the peak is strongly 
observable dependent and less strongly volume dependent. Without a monotonic behavior it 
is not possible to fit the r to a power law of the correlation length and therefore extract the 
dynamical critical exponent of the algorithm. This unexpected behavior is not an artifact of 
the algorithm used, since we observe it for multi-hit metropolis as well as with global HMC, 
but rather a feature of the physics of lattice QCD. Some author recently, in the framework 
of a different algorithm ||, observing the same behavior, conjectured some connection to 
the finite temperature phase transition. That interpretation is not fully consistent with 
previously well known results, about the position of that transition. From our analysis it 
is not possible today to offer support for that view since it is quite puzzling that the peak 
moves about according to the observables. Furthermore, even for the larger Wilson loops, 
the position of the peak is at values of (3 much lower than the value observed by the above 
authors. 
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Table 1: 

LHMC and Metropolis r for the 3x3 Wilson's loop. 
Table 2: 

t values for a variety of Wilson's loops. 

Table 3: 

Creutz's ratios. 



Table 4: 

Ratios of timings to cover one r. 
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